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Metastable states in Ising spin-glass models are investigated numerically by finding iterative so- 
lutions of mean-field equations for the local magnetizations rui . A number of iterative schemes are 
employed, and two different mean-field equations are studied: the Thouless- Anderson-Palmer (TAP) 
equations that are exact for the Sherrington-Kirkpatrick model, and the simpler 'naive-mean-field' 
(NMF) equations, in which the Onsager reaction term of the TAP equations is omitted and which 
are exact for the Wallace model. The free-energy landscapes that emerge are very different for 
the two systems. For the TAP equations, the numerical studies confirm the analytical results of 
Aspelmeier et al., which predict that TAP states consist of close pairs of minima and index-one 
(one unstable direction) saddle points, while for the NMF equations the corresponding free-energy 
landscape contains saddle points with large numbers of unstable directions. For the TAP equations 
the free energy difference between a minimum and its adjacent saddle point (the 'barrier height') 

scales as l/(/ — /o)^ where / is the free energy per spin of the solution and /o is the equilibrium 
free energy per spin. This means that for 'pure states' for which / — /o is of order 1 /N , where A'' is 

the number of spins in the system, the barriers between them scale as As, but between states for 
which f — fo is of order one, then the barriers are finite and also small so such metastable states 
will be of limited physical significance. For the NMF equations there are saddles of index K and we 
can demonstrate that their complexity Ek scales as a function of K/N. 

We have also employed an iterative technique with a free parameter that can be adjusted to bring 
the system of equations close to the 'edge of chaos'. Both for the TAP and NME equations it is 
possible with this approach to find metastable states whose free energy per spin is close to fo. As 
N, the number of spins is increased, it becomes harder and harder to find solutions near to the edge 
of chaos, but nevertheless the results which can be obtained are competitive with those achieved by 
more time-consuming computing methods and suggest that this method may be of general utility. 



PACS numbers: 75.10.Nr, 75.50.Lk, 05.70.Ln 
I. INTRODUCTION 



The study of metastable states in spin glasses began 
some 25 years ago with the calculations of Tanaka and 
Edwards "Q (TE) and of Bray and Moore (BM), deal- 
ing with metastable states at zero and non-zero temper- 
atures respectively, within the infinite-range Sherrington 
Kirkpatrick model. For the zero-temperature studies, 
metastable states were defined as states which are one- 
spin-flip stable, i.e. states for which flipping any one spin 
increases the energy, while for T > they were iden- 
tified with solutions of the TAP equations 0, which 
are exact for the SK model. While the TE calcula- 
tion is quite straightforward, the original BM calculation 
involved some technical challenges which were finessed 
in the first attempt, and have only recently been satis- 
factorily resolved, some 24 years after the original pa- 
per HM H- While the central result of BM, that the 
mean number (averaged over disorder configurations) of 
TAP solutions increases exponentially with the number 
of spins, {N,} - exp[NT,{T)], where E(r) is the 'com- 
plexity', was confirmed by the later work, new insights 
into the nature of the TAP solutions and the structure 



of the free-energy landscape were obtained. Specifically, 
it was found that all solutions correspond either to min- 
ima of the TAP free energy i^TAP, or to saddle points 
of index one, where the index of a saddle point is the 
number of negative eigenvalues of the Hessian matrix 

F-jAP / duiidrnj , so that minima have index zero. Fur- 
thermore, the minima and saddle points occur in close (in 
configuration space) pairs 5]. This structure has impli- 
cations for the dynamics in the SK model. In particular 
we will study the free energy difference between the sad- 
dle and its associated nearby minimum. This difference 
would seem intuitively to be related to the barrier which 
has to be overcome to escape from a pure state and we 
will find that it increases with TV like N^/^. However, for 
the vast majority of the TAP states — all those with free 
energy per spin larger than its equilibrium value — the 
barriers are finite. This means that such TAP states are 
probably of little physical significance for the dynamics 
of the SK model. 

We will contrast the TAP free-energy landscape with 
that of the naive mean-field (NMF) equations. We show 
numerically that the free-energy function i^NMF possesses 
saddle points with large indices K (up to a maximum 
proportional to TV), and we show that the correspond- 
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ing landscape is more rugged than the TAP landscape 
and tends to trap the 'iterative' dynamics at (or close 
to) a threshold free energy where the minima numeri- 
cally dominate the saddle points, much as in the p-spin 
spherical model 

This qualitative distinction between dynamics on the 
TAP and NMF free energy surfaces is interesting because 
both sets of equations become the same at T ~ 0. In par- 
ticular, they have the same ground states. This suggests 
a program for finding the ground-state energy by work- 
ing with the TAP equations at very low but non-zero 
temperature, and using the resulting states at starting 
configurations for a final quench at T = 0. In view of 
the preceding discussion, the same approach applied to 
the NMF equations would not be successful, since the 
system would get stuck at a threshold free energy, as for 
the p-spin spherical model. If one hopes to determine 
ground-state properties using T > equations, there- 
fore, it is important to choose equations with a favorable 
free-energy surface (like TAP) rather than an unfavor- 
able one (like NMF). In other words if one is trying to 
find low energy states of (say) the SK model using stan- 
dard algorithms such as gradient descent on a free energy 
surface, then one gets much lower energy solutions if one 
uses the exact free-energy surface than when one uses 
only an approximation to it. We suspect this point will 
have validity beyond the SK model. 

However, we have discovered that it is possible to ob- 
tain solutions of low free energy by a method which takes 
one entirely off the free energy surface. These states are 
accessed through a novel iterative algorithm containing 
one adjustable parameter. We show that there is a criti- 
cal value of this parameter separating runs which termi- 
nate from those which do not, and that the lowest free- 
energy states are accessed when the parameter is close to 
the critical value. As the parameter is increased towards 
the critical value, one sees many solutions which are limit 
cycles of a length which also increases as the critical value 
is approached and beyond the critical value only 'chaotic 
solutions' seem to exist. In the course of the iteration it 
is possible for the magnetizations to take unphysical val- 
ues such that I mi I > 1. We speculate that this approach 
may be a generic technique for finding low-cost states 
in hard optimization problems. Somehow being close to 
the 'edge of chaos' enables one to explore phase space 
and it is only the lowest states, which are surrounded by 
the highest barriers, which can trap the algorithm into 
a fixed point rather than a limit cycle. Possibly related 
ideas are in a paper by Boettcher and Frank j^. 

The outline of the paper is as follows. In section ^ 
we state the TAP and NMF equations, discuss the algo- 
rithmic methods we will adopt and demonstrate that for 
NMF most of the standard algorithms seem to terminate 
at a free energy which is relatively insensitive to the al- 
gorithm. In section llTTl we explain this result by carrying 
out an exhaustive search for turning points of the free en- 
ergy for small systems, showing that there is a threshold 
free energy below which minima numerically dominate 



saddle points. The connection between the structure of 
the free-energy landscape and the nature of the dynami- 
cal attractors is further explored in section HVI Sect ion M 
contains a detailed study of low free-energy TAP states 
obtained from studying the solutions using the 'edge of 
chaos' algorithm. Section IVIl concludes with a discussion 
and summary of the main results. 



II. SOLVING THE TAP AND NMF EQUATIONS 
BY ITERATION 



The TAP and NMF equations can be derived from the 
following generalized free-energy function: 



F = ~Y1 Jijm,mj - 7^(1 - if + 

(i.3) 
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where 7 = 1 corresponds to TAP and 7 = to NFM. 
In Eq. ^ the local magnetizations, rrii = (Si), lie in 
the range —l<mi<l,/3— 1/kT as usual, and q = 

The TAP and NFM equations themselves are derived 
from the equations dF/dnii = 0, i = 1,...,N, which 
give the turning points of F. They can be written in the 
form 



rrii = Gtim) = tanh [/3/ij - 7^^(1 - q)mt 



(2) 



where the local field hi — X^i^j JijUij. The second term 
in the argument of the tanh function is, for the TAP case 
7 = 1, the Onsager reaction term. The TAP equations 
are exact as — > 00 for the SK model. The naive mean- 
field equations are exact for the Wallace model, in which 
at each site i in the system there is a set of Ising spins 
Sia (a = 1, 2, . . . , A:), each of which interacts with the k 
spins at sites which are coupled to the i-th site and the 
limit A; — > 00 is taken Q. 

The set of spins rrii that satisfy 10) can be thought of 
as a fixed point fh* of the map 
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This map could be applied to all spins in parallel, to each 
spin separately and in sequence, each time to a random 
spin, or to a spin chosen for some particular reason. It 
turns out that for the NMF equations (7 = 0), applica- 
tion of this map leads always to a fixed point that coin- 
cides with a minimum in the free energy landscape, no 
matter how the map is applied. This fact one establishes 
by calculating, at the fixed point, the Hessian H which 
is the matrix of second derivatives of the free energy. 
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(which holds for all i,j if one takes Ju = 0) and checking 
that all its eigenvalues are positive, thus defining a min- 
imum. As we will see below, the free energy landscape 
also has turning points with a positive number K of neg- 
ative eigenvalues — we call these saddle points of index 
K. On the other hand the inclusion of Onsager reaction 
term leads to the map ^ almost never converging on a 
fixed point. 

Clearly one does not imagine that the trajectory de- 
scribed by iteration of the map ((SJ corresponds to any 
reasonable physical dynamics on the free energy land- 
scape, e.g., a zero-temperature dynamics in which one 
follows the path of steepest descent to a local (in gen- 
eral, metastable) minimum. Unfortunately, numerical al- 
gorithms that one might expect to approximate closely 
a physical dynamics (such as steepest descents or conju- 
gate gradients) turn out to be poorly adapted to the NMF 
free energy landscape. The reason for this is that the free 
energy minima tend to lie close to the sides of the hyper- 
cube — 1 < mi < 1 over which the free energy is defined. 
Naively, one might imagine that the divergence in the 
free energy gradient as the hypercube's boundary is ap- 
proached would prevent a descent algorithm from exiting 
the physical region. However, the weak and short-ranged 
logarithmic divergence is not resolved numerically. This 
can cause points to be reached that lie outside the hy- 
percube and for which the free energy is undefined, 
thereby creating an obvious numerical problem. We tried 
a number of modifications here to work around this, such 
as mapping the bounded hypercube to an infinite do- 
main, or inventing a ficticious free energy outside the 
hypercube, all to no avail. 

For this reason, we modified instead the iterative map 
(O in a way that we believe gives some insight into the 
free energy minima located by a physical dynamics. This 
map reads 



(fc+i) 



(fe) 
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in which a is a parameter that controls the speed of ap- 
proach to the next iterate. We describe the utility of 
this more general iterative scheme first in the context of 
the NMF equations (7 = 0) and then the TAP equations 
(7 = 1). 



iV 


Samples 


Random 


Cyclic 


Greedy 


Physical 


20 


2000 


-0.67196(3) 


-0.67151(3) 


-0.67285(3) 


-0.67211(3) 


30 


1000 


-0.68811(4) 


-0.68760(4) 


-0.68882(4) 


-0.68804(4) 


40 


1000 


-0.69900(4) 


-0.69862(4) 


-0.69973(4) 


-0.69892(4) 



TABLE L Average free energy of minima in the NMF free 
energy landscape reached using various iterative algorithms. 
For all system sizes N , /3 = 2 and different minima were 
located by restarting the algorithm for each sample from 5000 
different, random initial conditions. 



to lower the free energy, a fact that has been verified nu- 
merically. Although we do not claim that this dynamics 
is equivalent in any sense to a steepest-descent dynam- 
ics, the fact that the free energy is lowered in a controlled 
fashion leads us to view it as a 'physical' dynamics. 

In Tabled we display the free energy of the NMF min- 
ima located using four different algorithms, averaged over 
both a number of different realizations of the disorder 
(i.e., sets of the random variables , each of which we 
refer to as a sample) and over a number of different ini- 
tial conditions on the iteration. So that the algorithms 
can be properly compared for a given system size N , each 
was treated to the same set of samples. The system sizes 
used were quite small {N = 20,30,40) so that we may 
later compare with the exhaustive search for all turning 
points in the same samples (see Section ^Oj. The first 
three of the four algorithms are implementations of © 
applied on a spin- by-spin basis in different ways. In the 
random algorithm, the spin chosen for update is chosen at 
random, whereas in the cyclic algorithm, each is visited 
in order. Meanwhile, the spin chosen in the greedy algo- 
rithm is that for which the free energy will decrease the 
most under application of the map. Finally the physical 
algorithm uses Q applied simultaneously to all spins and 
with the control parameter a = 10~^. We notice from 
the Table that the differences in the free energy arising 
from a change of algorithm, whilst significant compared 
to the errors, are about two orders of magnitude smaller 
than those arising from an increase in the system size. 
Later, we will provide evidence that these changes are 
also small on the scale of free energies spanned by all the 
solutions of the NMF equations. 



A. NMF Equations 

If in Eq. the parameter a is small, one expects a 
smooth trajectory to be followed, and crucially one that 
is constrained to remain within the hypercube. This 
statement can be made more precise by noting that 
the displacement dnii = a(tanh/3/ii — rui) is bounded: 
\5mi\ < 2a and so can be made arbitrarily small by de- 
creasing a. Furthermore, one can show that the projec- 
tion of this displacement onto the gradient vector of the 
free energy is always negative. Thus, for sufficiently small 
a, the dynamics should take a small step in such a way as 



B. TAP Equations 

Solutions of the TAP equations, i.e. Q with 7=1, are 
much harder to find than those of the NMF equations. 
Part of the reason for this is that including the Onsager 
reaction term stabilizes the trivial (paramagnetic) solu- 
tion TO = 0. Indeed, when we employed simple-minded 
free energy minimization algorithms on the TAP free en- 
ergy landscape, we found almost always this trivial solu- 
tion. However, this solution is thermodynamically unsta- 
ble in the spin glass phase and so one must take steps 
to avoid it. One possibility, pursued elsewhere, is to use a 
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different free-energy function in the thermodynamically 
unstable region UJ. 

We succeeded in finding solutions of the TAP equa- 
tions using the iterative map (01 in the regime where a 
is greater than one. The reason why this choice is useful 
is that with a < 1 the iteration enters a limit cycle be- 
tween points that are unrelated to the turning points of 
the TAP free energy landscape. By amplifying the dis- 
placements between these points, it is possible to break 
out of the limit cycle. An added bonus is that paths to 
the trivial solution to = are also destabilized — this is 
because our observations suggest paths to this solution 
have an oscillatory character. 

Clearly, there is an optimum value for a for finding the 
largest number of solutions of the TAP equations. Too 
small, and the limit cycles are attractive; too large, and 
one may hop in great leaps about the hypercube, never 
to settle on a solution. As we shall see in Sec. the 
optimum a is in fact close to a point where the iterative 
dynamics become chaotic. As this point is approached, 
the dynamics do exit the hypercube, but there are no 
numerical problems associated with this since the right- 
hand side of (jSJ is defined even if Ito^I > 1. We only 
require that the trajectory followed by the iteration con- 
verges on a fixed point inside the hypercube in a finite 
time. As we shall also later explain in detail, when a 
is tuned to its optimum we can approach closely states 
with the thermodynamic equilibrium free energy. 



C. Preliminary Comparison of TAP and NMF 

We believe that, despite having the same T ^ limit 
as the TAP equations, solutions of the NMF equations 
will not take us to the lowest energy states. Unfortu- 
nately, since the equilibrium free energy of the NMF 
equations is not known exactly at finite temperature, we 
are unable to test whether it can be reached for this case. 
Instead, we probe the relationship between NMF free en- 
ergy minima and low-lying states by performing a zero- 
temperature quench from the former. Specifically, hav- 
ing found a solution to the NMF equations, i.e., ^ with 
7 = 0, at finite temperature /?, we then iterate the T = 
version of the map. That is, every spin is set to have unit 
magnitude but with the sign taken from the final state 
of the NMF iteration. Then, sites are visited randomly 
and flipped if they are not aligned with their local fields 
until a metastable state is reached in which all spins have 
become aligned with their fields. The energies of these 
states, averaged over different NMF turning points and 
different bond distributions, are shown in Fig. Ill CI for 
a range of system sizes N and inverse temperatures (3. 
A noticeable feature of the graph is that these energies 
have a nonmonotonic dependence on [3 and increase as 
the temperature is decreased past a certain point. This 
suggests that minima of the NMF free energy landscape 
(at least, those found by iterative methods) at low tem- 
peratures are far away from the low energy states. This 




FIG. 1: Energies reached from solutions of the NMF equations 
found by iteration after a quench to a T = state at different 
system size A'' and inverse temperature /3. 



one sees from the fact that the lowest energy states ob- 
tained by this process occur at an apparently well-defined 
finite temperature /9 ~ 1. 

A plausible explanation for the minimum is as fol- 
lows. The total number of NMF solutions has the form 
exp[iVI](r)], where S(T) is expected to be a monoton- 
ically decreasing function of T, vanishing at the critical 
temperature Tc — 2 (where there is only one solution, 
rui — for all i). Exactly this behavior is obtained 
for the TAP equations 0], except that = 1 for TAP. 
It follows that, for any given system size iV, there will 
be a range of temperatures near Tc where the condition 
NT,{T) ^ 1 will hold. In this regime there will be a 
unique solution with high probability. In fact there will 
be 0(1) solutions down to a temperature T*{N) defined 
by NT,{T*) = 1. For T > T* the free energy wiU be 
determined essentially exactly, as there is only one solu- 
tion. This state should provide a good starting point for 
the subsequent quench to T = 0. Our hypothesis, then, 
is that the minima in Figure ITl CI occur at T ~ T*{N). 
The relative insensitivity to N of the position of the min- 
imum in Figure III CI mav be accounted for by the strong 
T-dependence of S(T) in the regime T > Tc/2 (evident 
in the exact solution for TAP 

As we have mentioned above, solutions of the TAP 
equations by (suitably chosen) iteration at the 'edge of 
chaos' can as T ^ 0, yield states of the SK spin glass 
for which the energy per spin is close to the equilibrium 
value. We have also used an 'edge of chaos' algorithm 
for the T = version of the NMF equations and there 
too states of low free energy are obtained. However, the 
data presented in Fig. Ill CI strongly suggests that this is 
not the case for the NMF equations when treated by an 
iterative procedure which parallels gradient descents. We 
believe that this is a consequence of the shape of the free 
energy landscape, which we shall now discuss. (For the 
TAP equations gradient descents algorithms usually fail 
to give any solution other than the trivial solution with 
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all rrii = 0). 

III. THE SHAPE OF THE FREE ENERGY 
LANDSCAPE 

In this section we investigate numerically some proper- 
ties of the NMF free-energy landscape in order to contrast 
with a similar study of the TAP landscape presented in 
0. The strategy is to consider small systems and locate 
(as best one can) all the turning points in the landscape. 
We begin with a discussion of the numerical procedure. 




0.1 0.2 0.3 

A. Numerical method K/N 



In , turning points of the TAP free energy landscape 
were found using Broyden's method 12] for solving the 
system of nonlinear equations Q with 7 = 1. We found 
this algorithm also to be applicable to the case 7 = 0, 
albeit in a different implementation T3| . Although more 
sophisticated and efficient root-finding algorithms exist, 
we found that they tended to fail for much the same 
reasons as did steepest descent and conjugate gradient 
methods for locating minima (as discussed above). 

The numerical algorithm approaches a single root de- 
terministically from an initial condition, so in order to 
find different roots, a sequence of random initial con- 
ditions were tried. One hopes that, after a sufficiently 
large number of initial conditions, one will have located 
each of the turning points in the NMF landscape at least 
once. One test that reveals if any solutions are certainly 
outstanding is to calculate the Morse sum 

Sm = E Vl)^" (6) 
1=1 

where Ki is the saddle index (number of negative eigen- 
values of the Hessian) of the i*'' turning point. Topolog- 
ical considerations imply that Sm — 1 when the sum is 
over all turning points of the NMF free energy. Hence, 
if one finds Sm ^ 1, there remain solutions to be found. 
The converse, however, is not true, and so one must use 
a little trial-and-error to estimate the typical number of 
initial conditions after which no new solutions are found. 

Computational constraints imply a trade-off between 
temperature and system size. The larger the system size, 
the slower the search for a single solution and if one works 
at too high a temperature, this fact places too low a limit 
on the typical number of turning points per sample that 
can be located in a reasonable time. Conversely, at too 
low a temperature, the number of solutions per sample 
grows so fast with system size that one is unable to survey 
effectively a reasonable range of system sizes. We found 
an agreeable compromise occurred at /3 = 2 (the critical 
point in NMF is /3c = 1/2) which allowed the study of 
samples up to size N — 40. At this system size there were 
typically 600 turning points per sample, compared with 



FIG. 2: Complexity Yjk of the NMF equations with so- 
lutions separated according to the index of the saddle K. 
Near the pealc, the data collapse when K is scaled by TV" for 
iV = 20, 30, 40 at /? = 2. The dotted line corresponds to the 
index of the trivial solution in the thermodynamic limit. 



45 at = 20. Whilst the criterion Sm = 1 was satisfied 
for all but a handful of 2000 samples with A^ = 20 sites, 
it held only for about a third of the 1000 samples at A = 
40. The data presented in the following are obtained by 
amalgamating results from all samples, including those 
whose solution sets were certainly incomplete by virtue 
of Sm 7^ 1- Excluding these samples from the analysis 
changes our results very little. 



B. Results 

We first examine what classes of turning point arise in 
the free-energy landscape. In the TAP free energy land- 
scape, one finds only minima and index K = 1 saddles 
which occur in pairs (so that, along with the trivial solu- 
tion, the Morse sum is satisfied) From an analysis 
of all the minimum-saddle pairs for small systems, Cav- 
agna et al. ^ find that the free energy difference between 
a minimum and its corresponding saddle point seems to 
vanish in the thermodynamic limit. 

For the NMF equations we find minima and saddles 
of all indices up to some maximum that grows linearly 
with system size. This result is best illustrated using the 
saddle complexity E^e , defined as 

= (7) 

in which Nso\ns{K) is the number of solutions of the NMF 
equations that correspond to a saddle point of index K. 
Rescaling K with N one obtains Fig. |21 which shows a 
good collapse around the peak that occurs dX K k, O.IA. 

A little more insight is gained into the shape of the free 
energy landscape by plotting for a given saddle index K 
the complexity as a function of free energy — see Fig.|31for 



6 



the case iV = 40. Here one sees a set of roughly similarly- 
shaped curves with a well-defined maximum at a free 
energy that increases with K whilst the width of the 
curves decreases. This suggests that the highest-energy 
turning points are also the most unstable, and a closer 
inspection of the data reveals that the trivial solution 
seems always to have the greatest index. 

To determine the value of K/N of the trivial solution 
in the thermodynamic limit, we need to find the total 
density of negative eigenvalues of the Hessian matrix H, 
which from Q is, at m = (and 7 = for NMF), 

Hi] = -Jij + -j^Sij . (8) 

In this expression, J is a symmetric Gaussian random 
matrix, and one can appeal to the Wigner semicircle law 
to learn that, in the limit ^ cxd, its eigenvalues A have 
a distribution 

Pj(A)-^v/4-A2 (9) 

on the support — 2 < A < 2. Clearly, the diagonal part 
of H serves simply to shift the distribution to the right 
by 1//3 and so the number K of the N eigenvalues that 
are negative has the limit 

lini ^ = ^ f dXpj{X-l/P) (10) 

11 11 I T , , 

- 2-^"2^-4V'"^- ^''^ 

For (3 < 1/2 we see that the trivial solution has no 
negative eigenvalues, i.e., we are above the critical tem- 
perature for the spin glass phase. For the case (3 = 2 
studied numerically, the fraction of negative eigenvalues 
approaches K/N « 0.3425, which is plotted as a ver- 
tical dotted line in Fig. |21 We note that even for the 
small systems studied, the complexity approaches zero 
at this value. (Actually, in these finite systems, we ob- 
tained solutions with larger index but their complexity 
was negative — i.e. there was on average fewer than one 
such solution per sample — and hence irrelevant in the 
thermodynamic limit). 

More generally, one can look at the distribution of 
eigenvalues of the Hessian at other turning points. In 
Fig. 21 these are shown for the = 40 system, aver- 
aged over all turning points and also separated out into 
the average distributions at minima and saddles of index 
K = 1 and 2. Most noticeable is a dip in all the spectra 
around A = 0. We believe this is a real effect and not 
an artifact of the numerics: the plots for N = 20 — where 
we are confident that all the turning points in the free 
energy landscape have been isolated — show the same be- 
havior. We are however, currently unable to explain this 
repulsion of low-eigenvalue modes of the Hessian. 




FIG. 3: Complexity Ejf (/) of the NMF equations as a func- 
tion of free energy / for different saddle index K, system size 
iV = 40 and /? = 2.' 




FIG. 4: Eigenvalue distributions of the Hessian at minima, 
saddles of index 1 and 2 and all turning points. The sys- 
tem size = 40 and 13 = 2. Note that the vertical scale is 
logarithmic. 



IV. CONNECTION BETWEEN ITERATIVE 
SOLUTIONS AND THE STRUCTURE OF THE 
FREE-ENERGY LANDSCAPE 

We now move on to discuss how the structure of the 
free energy landscape infiuences (if at all) the attractors 
of the iterative maps described in the previous section. 
In the foregoing, we showed that for the NMF equations, 
the typical free energies of turning points of a particu- 
lar type increase with the saddle index K . In particular, 
there is a free energy at which, when approached from 
above, the number of minima exceed the number of sad- 
dle points. Recall, in fact, that the complexity describes 
the exponential increase of the number of solutions with 
system size, so that below this characteristic free energy 
the ratio of the number of saddle points to minima van- 
ishes in the thermodynamic limit. In Fig. |Slwe plot the 



N=20 





Free energy 



FIG. 5: Complexity T,K{f) of the NMF equations as a 
function of free energy / for three classes of turning point 
(minima, K = 1 saddles and all saddles) and system sizes 
A'' = 20, 30, 40, all at l3 — 2. For comparison, the mean free 
energies of states found by the various iterative algorithms 
are plotted. All the energies cluster together except for those 
obtained by the 'edge of chaos' algorithm i.e. Eq. @ with 
a = 1.9, which outperforms them by a considerable margin. 



mean free energies set out in Table for the various iter- 
ative algorithms alongside the corresponding complexity 
histograms for minima, K — 1 saddles and all saddles for 
the three system sizes N — 20, 30, 40. 

We first notice that the spread in free energies within 
the different iterative algorithms (indicating a slight bias 
of the algorithms towards different attractors) is small 
compared to the free-energy range spanned by a set of 
turning points (the minima, say). Noting that all three 
plots in Fig. [S] have the same horizontal scale, it is clear 
that as N increases, the mean free energy located by the 
iterative algorithms approaches the largest free energy 
at which the minima start to outnumber the saddles. Of 
course, since we have been unable to access larger system 
sizes this provides no proof that the these points will co- 
incide as N ^ oo; however, this data provides support 
for a 'landscape' picture of zero-temperature dynamics 
in which the free energy is lowered until such time as 
the minima dominate to such an extent that becoming 
trapped in a local minimum is an inevitability. Notice 
that in the 'saddles rule' type of hypothesis the dy- 
namics is supposed to be trapped at the value of the free 
energy at which the index density fc(/) = {K)/N be- 
comes non-zero, where the average is over all solutions 
whose free energy lies between / and f + df. We have 
plotted this quantity in Fig. El Clearly for the system 
sizes which we have been able to study the free-energy 
per spin of most algorithms is of order —0.69 (see Ta- 
ble m . This is much higher than the value of / at which 
k{f) appears to become non-zero. But the number —0.69 
does seem to fit better with the free energy at which in- 
dex 1 saddles begin to outnumber the index saddles i.e. 



FIG. 6: Index density of the NMF equations as a function of 
free energy f a,t (3 — 2 and system sizes — 20, 30, 40. 



the mimima, so perhaps a modified version of the 'saddles 
rule' hypothesis does apply to the NMF landscape. 

As mentioned above, the TAP free-energy landscape 
contains exponentially large (in N) numbers of minima 
and K — 1 saddles only. Not only that, but these turn- 
ing points come in pairs, i.e., a minimum that is close to 
a saddle in the sense that it is separated by a total free 
energy difference that vanishes approximately as N~^-'^^ 
in the limit — > oo @. What this means in particular 
is that for the TAP equations the complexity of minima 
So(/) and of K = 1 saddles Si(/) are equal for all free 
energies / (and all other complexities are zero). There- 
fore there is no 'crossing point' at which minima begin 
to outnumber saddles as there is for the NMF equations. 
Since the solutions we found are compatible with tracking 
down states of lower free energy, it would appear that the 
fact that there are equal numbers of saddles and minima 
over a range of free energy means that one does not get 
trapped in a particular minimum, but somehow manages 
to explore further the free energy landscape by passing 
over the nearby saddle-point. This latter statement is, 
however, as yet somewhat speculative and requires fur- 
ther investigation. 



V. DYNAMICAL PROPERTIES OF THE TAP 
ITERATION 

As mentioned previously, the TAP equations are usu- 
ally very difficult to solve numerically. However, using 
the modified map Eq. Q with a > 1 surprisingly leads 
to a vastly increased probability of finding solutions. In 
this section we explain in more detail why this is the 
case and how the ad hoc parameter a should be chosen 
for optimal results. First, however, we discuss the details 
of the implementation as there are a number of potential 
pitfalls. 
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A. Numerical iteration 

The iteration Eq. © is best applied to each spin sep- 
arately and in sequence to improve convergence. At low 
temperatures the converged result usually contains some 
rui which are equal to ±1 within numerical accuracy. 
While this might be considered merely an irrelevant nu- 
merical inaccuracy, it is a serious problem as soon as the 
Hessian, Eq. |0J, of the solution is required as it is clearly 
ill-defined if = ±1. Therefore we continue iteration 
of a transformed version of the TAP equations for the 
variables Xi = — sign(TOj) ln(l — mf ), 



(fc+i) (fc) , 
x) — X: + a 



+ 2sign(?i:,)ln 



1 + e 



-2\Hi\ 



xl 

(12) 



where Tii — tanh^^ Qi. This particular choice of trans- 
formed variables is adapted to the form of the Hes- 
sian, Eq. The diagonal elements of Hij can be ob- 
tained from \xi\ simply by exponentiation without fear 
of numerical overflow or loss of accuracy. The initial 
condition for the iteration of the transformed variables 
is obtained from the set of converged nii via xl"' = 

2Hi + 2sign(7ii) In "'""'"^ ^ ' applied to all spins in paral- 
lel. This is important since application in sequence would 
implicitly carry out more iteration steps and would not 
result in an accurate translation of the rrii into the trans- 
formed variables Xi, thus in practice often leading to exit 
from the basin of attraction. 

When convergence is reached using the transformed 
variables, the Hessian is calculated. Diagonalization of it 
is, however, hampered by its severe ill-conditionedness, 
particularly at low temperatures. The eigenvalues may 
vary over 100 orders of magnitude, owing to the ex- 
tremely large variation of the diagonal entries. Therefore 
standard diagonalization packages fail and one has to re- 
sort to, e.g., diagonalization routines from )15j modified 
for use with arbitrary precision packages such as |T^. 



B. Dynamical critical point 

The fact that iterations with a > 1 generally find many 
solutions of the TAP equations and iterations with a < 
1 basically never find solutions (other than the trivial 
one) suggests that a = ofc = 1 is a critical point of the 
iterative dynamics. Just as the iterative behavior of e.g. 
the logistic map is changed from convergence to a stable 
fixed point to eventually chaotic dynamics via its control 
parameter, we expect that we have stable fixed points 
for a > 1 and chaotic behavior for a < 1. We test this 
hypothesis by various means. First we analyze scaling 
properties which we expect to hold near the critical point. 
We then examine the lengths of limit cycles above and 
below the critical point. 
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FIG. 7: Finite-size scaling of the free energy as a function of 
the scaling variable x — N{a — ftc)^, with 7 = 1, at inverse 
temperature /3 = 20. 



C. Scaling properties near the critical point 

The choice of the parameter a is not arbitrary. A given 
value of a favors finding a particular subset of all solu- 
tions, i.e. the typical free energies found are correlated 
with a. Under the hypothesis of the existence of a criti- 
cal point at a = Qfc = 1 the free energies / obtained with 
a given a are hypothesized to obey a finite-size scaling 
form 



f-fo = N-^F{N{a-a,r), 



(13) 



where /o is the exact free energy per spin in the ther- 
modynamic limit, which is known exactly from methods 
described in [H [H to be /o « -0.76324 ... at /? = 20 
[l9i] . Fig. [7| shows the connection between a, the sys- 
tem size N and the average free energy / of the solu- 
tions found at inverse temperature /3 = 20 in a finite-size 
scaling plot with scaling variable x — N{a — ac)'^ . The 
best data collapse is obtained for an exponent 7 ~ 1, but 
with too much scatter to provide an error bar. The finite- 
size corrections are still very large, as can be seen from 
the deviations at large arguments. The scaling function 
F{x) seems to behave as x"^ for large x with 77 « 1. This 
would imply that, away from the critical point, the excess 
free energy per spin approaches an algorithm-dependent 
(here a-dependent) value, in accordance with a result of 
Newman and Stein |23| which states that, in the thermo- 
dynamic limit, any given algorithm for solving the TAP 
equations will find solutions of a given free energy per 
spin, characteristic of the algorithm. In our case any 
fixed choice of a defines an algorithm and the scaling 
form with r; = 1 then gives f — fo + const. (q — ckc) for 
N — + 00, as required. In the opposite limit, the data 
suggest -F(O) = const., i.e. for a ac, the excess free 
energy per spin vanishes. 

As a spin-off of these results we therefore deduce that 
using the modified map Eq. we can, in principle, find 
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FIG. 8: The data from Tab. HH plotted for lengths 1, 2, and 



FIG. 9: Success rates (probability of finding a solution) as a 
function of the scaling variable N{a — Qc). 



states with free energy per spin arbitrarily close to /o by 
choosing the scaling variable x — N{a — Uc)^ as small as 
possible. The drawback is that the probability of finding 
any solutions at all goes to zero as a —> ac. This will be 
illustrated in the following two subsections. 



D. Cycle lengths 

As a is decreased, limit cycles in the dynamics occur 
with increasing frequency. In order to further test the hy- 
pothesis that a = 1 is a critical point we have calculated 
the lengths of limit cycles as a function of a. Table Hll 
and Fig. IHl show the frequency with which limit cycles of 
a given length were found for 10 samples of an = 800 
system at /3 = 20 using 25 different random starting po- 
sitions each. We checked for limit cycles up to length 20 
beyond which they were considered to be infinitely long. 
There is a marked change from a = 1.1, where 65% of the 
starting positions lead to convergence (i.e. a limit cycle 
of length 1), to a = 1, where 100% of the starting posi- 
tions did not run into a limit cycle (of length less than 
20) at all. There is another very distinct change between 
a — 0.99 and a = 0.95 from infinite limit cycles to cycles 
of length 2. Closer inspection revealed that these latter 
cycles are of the form rrii — > —rrii rrii. 

The picture which emerges from this data is that for 
a > 1 we find a distribution of cycle lengths which in- 
cludes a finite fraction of cycles with length 1, i.e. con- 
vergence. When a is tuned towards 1, more and more 
weight is taken from the finite limit cycles and trans- 
ferred to the infinite one, until at a = 1 all cycles have 
infinite length. At the same time, decreasing a towards 
the value one leads to states with ever lower free ener- 
gies. The optimum results for the free energy (giving the 
equilibrium free energy per spin in the thermodynamic 
limit) seemingly occur when the system is 'at the edge of 
chaos'. Similar ideas have put forward in a recent paper 
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FIG. 10: Success rates as a function of ^^^{a — Uc). 

by Boettcher and Frank who refer to 'optimizing at 
the ergodic edge.' 



E. Success rates 

The probability of finding cycles of length 1, i.e. con- 
vergence, is not only a function of a but, naturally, also 
of the system size A^ . While the scaling of the free energy 
in Fig.[71suggests that by tuning a towards 1 solutions of 
arbitrarily low free energy can in principle be found, we 
will now show that in practice this becomes increasingly 
difficult for large system sizes. Fig. |51 shows the success 
rates (probability of convergence to a solution from a 
random starting position) for various system sizes. The 
larger the system is, the more attempts are needed to find 
a solution. This statement can be made more precise by 
plotting the success rates against x' = N~^/^{a — ac), 
which is done in Fig. ^| It shows that, at least for sys- 
tem sizes > 800, the success rates fall on a master curve 
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TABLE 11: Frequency of appearance of limit cycles of a given length. The system size is = 800 at temperature /3 — 20. 
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FIG. 11: Finite-size scaling of the barrier height as a function 
of the scaling variable x — N{a ~ ac) at inverse temperature 
(3 = 20. 



FIG. 12: The (total) barrier height as a function of free energy 
(per spin). Each data point represents the average over all 
solutions found with a particular value of a. 



which goes to zero very quickly (faster than a power, and 
consistent with a Lifshitz-hke behavior ^ g-const.xa; 
with C « 2) for small arguments. 

So while the appropriate scaling variable is a: = iV(Q; — 
ac) for the free energies and barrier heights (see below), 
it is x' — N^^/'^x for the success rates. Fixing x and thus 
a preferred range of free energies drastically reduces the 
success rates for large systems. Assuming a Lifshitz tail 
as suggested above, the success rates are proportional to 

g-const.xAr''</^ 



F. Barrier heights 

Since minima and saddles of index 1 always occur in 
pairs in the TAP landscape 0, the barrier height of a 
solution can be defined as the free energy difference be- 
tween a minimum and its corresponding saddle. Fig. 1111 
shows a scaling plot of the barrier heights using the same 
scaling variable as in Fig. [7| The data collapse is not 
very good, but is consistent with a barrier height expo- 
nent ip = 1/3, which is the expected scaling for the SK 



model |2l| . This suggests that the barriers between very 
low-lying mimima and their corresponding saddles, with 
free energies per spin equal to the equilibrium value, scale 
as N'^/^ for large N ^ but with a rather small coefficient 
of order 10~^. For larger values of the scaling variable, 
where the free energy per spin lies above the equilibrium 
value (see Fig.O, the barriers decrease. The straight line 
in Fig. ^2 has a slope —1/3, so if the data were to fol- 
low this line it would indicate that the barriers become 
A^-independent at higher free energies. Recall that the 
barriers averaged over all TAP states were found to ac- 
tually decrease with A^, roughly as N~^-'^^ 0. It should 
be noted, however, that the states of a given free energy 
reached by our particular algorithm almost certainly do 
not sample all states of that free energy uniformly. 

The scaling behavior of the excess free energy per spin 
and of the barrier heights for large argument x — N{a — 
Q!c), namely f ~ fo ^ x and N~^/^B ~ a:^^/"^, indicates 
that (at least for large x) the barrier height is directly 
proportional to (/ — fo)~^^^ with no dependence on N 
or a — ac- Fig. 1121 demonstrates that surprisingly this 
is not only true for large x but for all values of A^ and 
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FIG. 13: Fini te size scaling plot of the average smallest eigen- 
value Xmin{q) of the Hessian at a minimum versus 1 — q. The 
eigenvalues scale with N~^^^ . The data shown is for a such 
that the scaling variable is a constant, namely N{a—ac) — 16. 



a — ac explored. This is an important result since the 
dependence on the unphysical a has been eliminated and 
the relation B ^ (f — /o)~^^^ is a property of physical 
quantities alone. 



G. Eigenvalues of the Hessian 

According to our picture of the TAP landscape de- 
scribed in [5j, the Hessian at a minimum should, in the 
thermodynamic limit, have one null eigenvalue and a 
band of eigenvalues starting at a value strictly larger 
than zero. In order to confirm this, we have studied the 
eigenvalues of the Hessian at the solutions found by our 
algorithm. Since the arbitrary-precision diagonalization 
is forbidding for large system sizes, we are restricted to 
N < 800. In the range of accessible system sizes we 
find a huge variation in the smallest eigenvalue Amin (see 
Fig. I13|l which is at first sight contradicting our land- 
scape picture. However, the eigenvalues Amin show a 
strong correlation with the corresponding value of q at 
the minimum. This is understandable from the form of 
the Hessian, Eq. since the smallest eigenvalue will be 
of the same order of magnitude as the smallest diagonal 
entry of Hij, which is proportional to 1/(1 — m^j^j) {rrinun 
being the m* closest to 0). Therefore, for q very close to 
1, Amin ~ 1/(1 — ™min) ~ 1/(1 ~ l) ■ Fluctuations of q 
towards 1 will therefore cause great fluctuations in Amin- 
In the thermodynamic limit, q should be self-averaging, 
therefore the variations in Amin become small and our 
observation merely indicates the enormous finite-size cor- 
rections still present at our accessible system sizes. 

The correlation between Amin and q and the strong 
variation of Amin over several orders of magnitude indi- 
cate that a simple average of Amin might not be very 
informative. Instead we consider the average Amin('z) of 



FIG. 14: Overlap between the eigenvector emin of the Hessian 
belonging to the smallest eigenvalue and the normalized dif- 
ference vector Am between a minimum and the corresponding 
saddle as a function of Ag — Igsaddio — gmin|. The data are for 
iV = 800 at /3 = 20 and a = 1.02. 

the eigenvalues as a function of q. This quantity is plot- 
ted in Fig. ^] in a finite size scaling plot which shows 
good data collapse, in particular for larger TV. There are 
several points to notice in this figure. First, the eigenval- 
ues scale with N~^/^ but can still be very large (of the 
order 1000 or more at the very low temperature, (3 — 20, 
at which we are working). Second, for larger N the range 
of g- values shrinks. This must be the case as in the ther- 
modynamic limit there should be only one value of q. 
Third, the range of Amin also shrinks and Amin (where 
the average now includes all values of q) goes to zero 
as N~^/^ (note that it would have been very difficult to 
obtain this result from a simple g-independent average 
of Amin)- This is in support of the picture of the TAP 
free energy landscape outlined previously which as- 
serts that minima and saddles occur in pairs which move 
closer together as the system size grows. The smallest 
eigenvalues of the Hessian therefore have to go to zero 
with increasing N. 

In order to provide further support for this picture, we 
have plotted in Fig. ^] the overlap between the eigen- 
vector Cmin belonging to the smallest eigenvalue and the 
normalized difference vector (in m-space) between the 
minimum and its corresponding saddle Am, as function 
of the difference Aq = | (/saddle — 9min|- If the two vec- 
tors were uncorrelated, the distribution of their overlaps 
would be sharply peaked around zero. The figure shows, 
however, that they are strongly correlated and that the 
direction of the smallest eigenvector coincides with the di- 
rection in which the corresponding saddle is found. The 
correlation is particularly strong for small differences of 
q between saddle and minimum since these pairs are very 
close together. 

The second smallest eigenvalue A2 is expected to be, in 
the thermodynamic limit, strictly larger than 0. In order 
to test this, we have plotted in Fig. E| the difference 
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FIG. 15: Plot of the difference between the second and the 
lowest eigenvalue as a function of A^^'^^Amin- 
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FIG. 16: Finite-size scaling of the NMF free-energy at T = 
as a function of the scaling variable x — N{2 — a). 



between A2 and Amin as a function of A^^/^Amin. In the 
thermodynamic hmit there should be only a horizontal 
line of finite width in this plot. Clearly we are still very 
far away from this limit, even for N — 800, the variation 
in A2 — Aniin being of order 10 or more. However, it is also 
clear that the data becomes more and more concentrated 
for larger N. Observe also the emergence of a horizontal 
line for N = 800. 



VI. DISCUSSION AND SUMMARY 

In this work we have compared and contrasted the 
free-energy landscapes of two closely related models, de- 
scribed by the TAP equations and the NMF equations 
respectively. The TAP free-energy surface has a sim- 
ple structure with turning points of only two types — 
minima, and saddle points of index 1. Furthermore, the 
minima and saddle-points occur in pairs, each containing 



a minimum and a saddle point that are close neighbors 
in configuration space. By using a novel algorithm con- 
trolled by a single parameter a, we are able to access TAP 
states whose free energy per spin is equal to the equilib- 
rium free energy per spin at the given temperature, with 
an excess total free energy of order 10 (see Fig.[7J|. This is 
achieved by choosing the parameter a close to a critical 
value, ttc, beyond which the iterative dynamics becomes 
chaotic. The dependence of the excess total free energy 
AF on N and a is broadly consistent with a scaling form 
AF — F[N{a ~ ac)'^], valid in the limit N — > c», a ac 
with N{a — ac)'^ fixed, and 7 « 1. The scaling function 
F{x) should behave as F{x) ~ x for large x, so that AF 
becomes of order N away from the critical point, while 
F(0) is a non-zero constant. 

We have also studied the free-energy barriers, given 
by the free-energy difference between a minimum and its 
corresponding saddle point. In the low free-energy regime 
where N{a—ac)'^ <^ 1, the barriers seem to grow as N-^^^, 
albeit with a small coefficient, while for N{a — ac)^ ^ 1 
we expect the barriers to be of order unity or smaller (see 
Fig. Ill|l . This leads to a picture where only the lowest 
free energy states (with free energies per spin equal to the 
equilibrium one) are separated by barriers which increase 
in some manner with N. 

In contrast to the TAP free-energy function, the NMF 
free energy has saddle points with index K taking val- 
ues up to order TV. The separate complexities £;<■(/) 
have been computed numerically (see Fig. |3J). The free 
energy at which the minima start to outnumber (expo- 
nentially in N) the saddle points seems to play the role 
of a threshold free energy which acts as lower bound on 
the free energy that is dynamically accessible (Fig.jSJ. 

We conclude with with some general remarks concern- 
ing optimization problems. If one is faced with a combi- 
natorial optimization problem, such as finding the ground 
state of a spin glass, one approach is to write down TAP- 
like equations and solve them using an iterative scheme 
which is related to a descent procedure on the free-energy 
landscape. As the temperature is lowered the states ob- 
tained approach those of a zero-temperature problem. In 
the spin-glass context, it does not seem to matter, a pri- 
ori, whether the finite-temperature generalization of the 
problem is through the TAP equations or the NMF equa- 
tions, since they become identical at T = 0. We believe, 
however, that the very different structures of the free- 
energy landscapes in the two models will make a very big 
difference in practice. While the NMF free-energy surface 
is 'robust', with its threshold free energy and, presum- 
ably, large free energy barriers between states, the TAP 
free energy surface may be termed 'fragile'. With any 
such problem therefore, it will pay to construct finite- 
temperature equations corresponding to a free energy 
with a fragile structure. This process might be called 
'sculpting the free-energy landscape'. 

Finally we would like to point out that using algo- 
rithms which take one off the free-energy surface into 
physically unreachable regions before reaching the final 
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physical fixed point may also have great utility, as illus- 
trated by our 'edge of chaos' algorithm. The success of 
this type of algorithm seems to be partly due to the fact 
that by allowing unphysical values of the one can tun- 
nel through barriers and partly because tuning the ad- 
justable parameter a towards the edge of chaos yields the 
solutions of low free energy which tend to be surrounded 
by the highest barriers which can trap the algorithm into 
a fixed point solution. We mainly used the 'edge of chaos' 
algorithm in connection with solving the TAP equations 
but it seems to work well also for the NMF equations. In 
Fig. [3 we have included the results obtained by setting 
a — 1.9 which is close to the 'edge of chaos' but with- 
out attempting to tune it to its optimum value for each 
value of N . (For the NMF equations the edge of chaos 
seems to be close to 2 and the values of a which give sta- 
ble solutions are less than the critical value). It is clear 
that it outperforms the other iterative algorithms, which 
approximate to descents on the free-energy surface, by 
a considerable margin. At T = the algorithm is espe- 
cially simple and Fig. 1161 shows that it seems to work as 
efficiently as the 'edge of chaos' algorithm applied to the 
TAP equations at low temperatures. However, it does 
suffer from similar problems as the 'edge of chaos' algo- 
rithm when applied to the TAP equations: we noticed 
that the number of iterations required to achieve con- 
vergence to a fixed point increased dramatically as the 



system size N increases, or as a approaches its critical 
value. 

In the context of mean-field problems such as finding 
the groundstate of the SK model, the 'edge of chaos' 
algorithm works well. We have also tested it on the one- 
dimensional spin glass at T = 0, (which can of course 
be solved exactly, enabling one to easily judge the ac- 
curacy of any proposed treatment). Although it again 
outperformed simple descent algorithms it was hard to 
get solutions whose energies per spin were equal to the 
equilibrium energy per spin. We attribute the difference 
to the nature of the free-energy landscape. In one di- 
mension the barriers around the states of the lowest en- 
ergy are little different in height to those surrounding 
metastable states of much higher free energy. 
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